function [ll, grad] = T1IRAbsParameters( par, x, echoes, sigma, scales, TR)

% computes the likelihood of the rice distributed inversion recovery data
% echoes = inversion times [n x 1] vector
% par = the 3 element parametervector [A, B, R1]
%
% The datamodel is given by
% A*(1-B*exp(TI*R1) + exp(TR*R1) )
% M = A/scales(1) - (B/scales(2))*1000*exp( - TI * (R1/scales(3)) ) + rice distributed noise with sigma = 1.
% dM/dPar = {1, -exp(-TI * R1), TI*B*exp(-TI * R1)}
% B_0 = R2, echoes = TEs
%
% Created by Henk Smit, EMC, 01-2011 based on the work by Dirk Poot, University of Antwerp, 13-8-2007

numTE = size(echoes,1);

if (size(par,1) ~= 3 || size(x,1)~=numTE)% + (nargin<=3)  HENK
    error('incorrect par or TE input');
end;

par = par.*scales;
A = par(1);
B = par(2);
R1 = par(3);
t = echoes;

% E1 = exp(min(650, -t*R1));
% S1 = A - B*E1;
% S = abs(S1);

% %mckinzey with TR
E1 = exp(min(650,-t*R1));
E2 = exp(min(650, -(TR)*R1));
S1 = A*(1-B*(E1)+E2);
S = abs(S1);

if nargout>1   
    sgn = sign(S1);
    
%     dSdR1 = (t*B) .*E1 .* sgn;
%     dSdA  = sgn;
%     dSdB  = -E1.*sgn;
%     dSdpar = [dSdA dSdB dSdR1];
    
%     %mckinzey with TR
    dSdR1 = ((t*A*B) .*E1 - A.*(TR).*E2) .* sgn;
    dSdA  = (1 - B.*E1 + E2) .* sgn;
    dSdB  = -A.*E1.*sgn;
    dSdpar = [dSdA dSdB dSdR1];
    
    [lrpdf, ricegrad] = logricepdf(x, S, sigma,logical([0 1 nargin<=3]));
    
    grad = -sum(bsxfun(@times,ricegrad,dSdpar))';
    grad = grad.*scales;
else
    [lrpdf] = logricepdf(x, S, sigma );
end;
ll = - sum( lrpdf(:) );
if ~isfinite(ll)
    ll = inf;
    grad = zeros(size(grad));
end;